import unittest
import numpy as N
import ftreactormodel
import ratelaw
import pdes
import pellet

class ModelTestRunner(unittest.TestCase):
    def testftmodel(self):
        bed = [range(25), range(25)]
        feed = [0.1, 0.2, 0.3, 0.4, 0.5, 123, 234]
        ratelaw = [1, 2]
        material = [1, 2]
        solveropt = [1, 2]
        ftreactormodel.ftreactormodel(bed, feed, ratelaw, material, solveropt)
        
    def testlhratelaw(self):
        params = range(22)
        r = ratelaw.LH2RateLaw(params)
        r.calc_rate([1, 2], 300)
        
    def testbaseratelaw(self):
        params = range(4)
        r = ratelaw.RateLaw('DETAIL', [1.1, 2.1], 1)     
        self.assertRaises(NotImplementedError, r.calc_rate, [], [])

class PDETestRunner(unittest.TestCase):
    def testtrisys(self):
        a = [3., 2., 1.]
        b = [1., 7., -2., 5.]  # constant vector
        c = [4., 1., 4.]
        d = [1., 4., 3., 3.]
        r = pdes.trisys(a, d, c, b)
        print(r)

class PelletTestRunner(unittest.TestCase):
    def _testpelletlh2(self):
        """Test L-H -2 ratelaw"""
        params = [549,86000,.11,89000,.1,0,.5,.75,.5,.25,0,2,1,0,0,1.25,1.25,.25,0,1,-160000,-210000]
        ra = ratelaw.get_ratelaw('L-H-2', params)
        radius = 0.5/10/2
        shellNum = 20
        p = pellet.Pellet(ra, radius, 0.98, N.array([1e-6, 1.2e-5, 1e-6, 1e-6]), 0.01, 0.01)
        propArray = N.array([0.315222466020354, 0.627011329024444, 0.0013, 0.02])
        raArray, effCO, sels, deltaH, profileArray = p.calc(True, 20*1.013e+5, propArray, 218+273.15, shellNum, 1)
        print 'raArray:'
        print raArray
        print 'sels:'
        print sels
        print 'effCO:'
        print effCO
        
        import pylab
        x = N.linspace(0.0, radius, shellNum)
        pylab.plot(x, profileArray[0], 'b-', x, profileArray[1], 'b-',
                   x, profileArray[2], 'b-', x, profileArray[3], 'r-')
        pylab.title('$LH2Concentration$')
        pylab.show()
    
    def _testpelletdetail(self):
        """Test DETAIL rate law"""
        carbonNum = 50
        params = [carbonNum,.00000081,1152,92.89e+3,113.76,87.01e+3,50400,111.4e+3,0.0,45.08e+3,.0000275,.00113,.0181,4.68,22.6,-160000,-210000]
        ra = ratelaw.get_ratelaw('DETAIL', params)
        radius = 1.0/10/2
        shellNum = 20
        p = pellet.Pellet(ra, radius, 0.98, N.r_[1e-6, 1.2e-5, 1e-6, 1e-6, 1e-6, 1e-6 * N.ones(2 * (carbonNum - 1))], 0.01, 0.01)
        propArray = N.r_[0.315222466020354, 0.627011329024444, 0.0013, 0.02, 0.0, N.zeros(2 * (carbonNum - 1))]
        raArray, effCO, sels, deltaH, profileArray = p.calc(True, 20*1.013e+5, propArray, 218+273.15, shellNum, carbonNum)
        print 'raArray:'
        print raArray
        print 'sels:'
        print sels
        print 'effCO:'
        print effCO
        
        import pylab
        x = N.linspace(0.0, radius, shellNum)
        pylab.plot(x, profileArray[0], 'b-', x, profileArray[1], 'b-',
                   x, profileArray[2], 'b-', x, profileArray[3], 'r-')
        pylab.title('$DETAILConcentration$')
        pylab.show()
        pylab.save('pdetail.png')
        
    def test_ratio_variation(self):
        """Variation of H2/CO ratio along the pellet dimension
        H2/CO = 1.99
        T = 218 + 273.15 K
        Rp = 1.0 mm
        Pressure = 15, 20, 25 atm
        """
        carbonNum = 25
        params = [carbonNum,.00000081,1152,92.89e+3,113.76,87.01e+3,50400,111.4e+3,0.0,45.08e+3,.0000275,.00113,.0181,4.68,22.6,-160000,-210000]
        ra = ratelaw.get_ratelaw('DETAIL', params)
        radius = 2.0/10/2
        shellNum = 20
        pressure = [15, 20, 25]
        data = N.zeros((len(pressure), shellNum))
        p = pellet.Pellet(ra, radius, 0.98, N.r_[1e-6, 1.2e-5, 1e-6, 1e-6, 1e-6, 1e-6 * N.ones(2 * (carbonNum - 1))], 0.01, 0.01)
        propArray = N.r_[0.315222466020354, 0.627011329024444, 0.0013, 0.02, 0.0, N.zeros(2 * (carbonNum - 1))]
        
        for n in range(0, len(pressure)):
            raArray, effCO, sels, deltaH, profileArray = p.calc(True, pressure[n]*1.0e+5, propArray, 218+273.15, shellNum, carbonNum)
            data[n] = profileArray[1] / profileArray[0]
        
        import pylab
        x = N.linspace(0.0, radius, shellNum)
        pylab.semilogy(x, data[0], 'bs-', x, data[1], 'bo-', x, data[2], 'b^-')
        pylab.xlabel('Pellet Dimension')
        pylab.ylabel('H2/CO fugacity ratio')
        #pylab.legend(str(pressure[0]) + 'bar', str(pressure[1]) + 'bar', str(pressure[2]) + 'bar')
        pylab.show()
        
            
        
        